#include "texArrayHandle.h"
#include "kernel.h"
#include "molecule.h"
#include "parse.h"
#include <map>
#include <utility>
#include <bitset>
#include <math.h>
#include <stdio.h>

using std::map;
using std::pair;
using std::bitset;

#define MAX_ATM 1000
#define MAX_TOR 50

extern const int IDX_E=200;
extern const int IDX_D=300;
extern const int IDX_VDW=500;
extern const int IDX_SOL=501;
extern const int IDX_EPS=502;
extern const int IDX_R_EPS=503;

texture <float,3,cudaReadModeElementType> tex_H;	//0
texture <float,3,cudaReadModeElementType> tex_HD;	//1
texture <float,3,cudaReadModeElementType> tex_HS;	//2
texture <float,3,cudaReadModeElementType> tex_C;	//3
texture <float,3,cudaReadModeElementType> tex_A;	//4
texture <float,3,cudaReadModeElementType> tex_N;	//5
texture <float,3,cudaReadModeElementType> tex_NA;	//6
texture <float,3,cudaReadModeElementType> tex_NS;	//7
texture <float,3,cudaReadModeElementType> tex_OA;	//8
texture <float,3,cudaReadModeElementType> tex_OS;	//9
texture <float,3,cudaReadModeElementType> tex_F;	//10
texture <float,3,cudaReadModeElementType> tex_Mg;	//11
texture <float,3,cudaReadModeElementType> tex_P;	//12
texture <float,3,cudaReadModeElementType> tex_SA;	//13
texture <float,3,cudaReadModeElementType> tex_S;	//14
texture <float,3,cudaReadModeElementType> tex_Cl;	//15
texture <float,3,cudaReadModeElementType> tex_Ca;	//16
texture <float,3,cudaReadModeElementType> tex_Mn;	//17
texture <float,3,cudaReadModeElementType> tex_Fe;	//18
texture <float,3,cudaReadModeElementType> tex_Zn;	//19
texture <float,3,cudaReadModeElementType> tex_Br;	//20
texture <float,3,cudaReadModeElementType> tex_I;	//21
texture <float,3,cudaReadModeElementType> tex_Z;	//22
texture <float,3,cudaReadModeElementType> tex_G;	//23
texture <float,3,cudaReadModeElementType> tex_GA;	//24
texture <float,3,cudaReadModeElementType> tex_J;	//25
texture <float,3,cudaReadModeElementType> tex_Q;	//26

texture <float,3 ,cudaReadModeElementType> tex_e;	//200
texture <float,3 ,cudaReadModeElementType> tex_d;	//300

texture <float,3 ,cudaReadModeElementType> tex;

texture <float,3, cudaReadModeElementType> vdw_tables;
texture <float,1, cudaReadModeElementType> sol_ruler;
texture <float,1, cudaReadModeElementType> eps_ruler;
texture <float,1, cudaReadModeElementType> r_eps_ruler;

//__device__ __constant__ float sol_ruler[1000];
//__device__ __constant__ float r_eps_ruler[1000];

__device__ __constant__ float coord_x[MAX_ATM];
__device__ __constant__ float coord_y[MAX_ATM];
__device__ __constant__ float coord_z[MAX_ATM];
__device__ __constant__ float elec_q[MAX_ATM];
__device__ __constant__ int type[MAX_ATM];
__device__ __constant__ bool ignore[MAX_ATM];
__device__ __constant__ int atomCount;
__device__ __constant__ int flexCount;
__device__ __constant__ int ligandCount;
__device__ __constant__ int torCount;
__device__ __constant__ int torsion[MAX_TOR][2];
__device__ __constant__ unsigned long move[MAX_ATM];
__device__ __constant__ float axn_x[MAX_TOR];
__device__ __constant__ float axn_y[MAX_TOR];
__device__ __constant__ float axn_z[MAX_TOR];

__device__ __constant__ float VOL[27]; // hard-coded here, 27 atom types in total;
__device__ __constant__ float SOLPAR[27]; // hard-coded here, 27 atom types in total;


__device__ __constant__ float min_coord[3];
__device__ __constant__ float spacing;

//extern const int IDX_E=200;
//extern const int IDX_D=300;

extern "C" void mapSet(float &step, float *min)
{
	cudaMemcpyToSymbol(spacing,&step,sizeof(float));
	cudaMemcpyToSymbol(min_coord,min,sizeof(float)*3);
}

extern void init_texture(tex3DHandle<float> &handle, int index,bool filtered)
{
	const textureReference * texPtr;
	if      (index==IDX_E) cudaGetTextureReference(&texPtr, "tex_e");
	else if (index==IDX_D) cudaGetTextureReference(&texPtr, "tex_d");
	else if (index==IDX_VDW) cudaGetTextureReference(&texPtr, "vdw_tables");
	else if ((index>=0)&&(index<=26))cudaGetTextureReference(&texPtr, texName[index]);
	textureReference * texRefPtr=const_cast<textureReference *> (texPtr);
	texRefPtr->normalized=false;
	if (filtered)
		texRefPtr->filterMode=cudaFilterModeLinear;
	else
		texRefPtr->filterMode=cudaFilterModePoint;

	texRefPtr->addressMode[0]=cudaAddressModeClamp;
	texRefPtr->addressMode[1]=cudaAddressModeClamp;
	texRefPtr->addressMode[2]=cudaAddressModeClamp;

	cudaBindTextureToArray(texPtr,handle.getPtr(),handle.getChannel());
}

extern void init_texture(tex1DHandle<float> &handle, int index)
{
	const textureReference * texPtr;
	if (index==IDX_SOL) cudaGetTextureReference(&texPtr, "sol_ruler");
	else if (index==IDX_EPS) cudaGetTextureReference(&texPtr, "eps_ruler");
	else if (index==IDX_R_EPS) cudaGetTextureReference(&texPtr, "r_eps_ruler");
	textureReference * texRefPtr=const_cast<textureReference *> (texPtr);
	texRefPtr->normalized=false;
	texRefPtr->filterMode=cudaFilterModePoint;
	texRefPtr->addressMode[0]=cudaAddressModeClamp;

	size_t offset=0;

	cudaBindTexture(&offset,texPtr,(void*)handle.getPtr(),handle.getChannel(),1000*sizeof(float));
}

extern "C" void ligandCopy(ligand &lig)
{
//const float vol[]={0.0000,0.0000,0.0000,33.5103,33.5103,22.4493,22.4493,22.4493,17.1573,17.1573,15.4480,1.5600,38.7924,33.5103,33.5103,35.8235,2.7700,2.1400,1.8400,1.7000,42.5661,55.0585,33.5103,33.5103,33.5103,33.5103,33.5103};
//const float solpar[]={0.00051,0.00051,0.00051,-0.00143,-0.00052,-0.00162,-0.00162,-0.00162,-0.00251,-0.00251,-0.00110,-0.00110,-0.00110,-0.00214,-0.00214,-0.00110,-0.00110,-0.00110,-0.00110,-0.00110,-0.00110,-0.00110,-0.00143,-0.00143,-0.00052,-0.00143,-0.00143};
	cudaMemcpyToSymbol(VOL,&vol,sizeof(float)*27);
	cudaMemcpyToSymbol(SOLPAR,&solpar,sizeof(float)*27);

	int n=lig.atomCount;
	int t=n-lig.flexCount;
	cudaMemcpyToSymbol(atomCount,&lig.atomCount,sizeof(int));
	cudaMemcpyToSymbol(flexCount,&lig.flexCount,sizeof(int));
	cudaMemcpyToSymbol(ligandCount,&t,sizeof(int));
	cudaMemcpyToSymbol(coord_x,&lig.coord_x[0],sizeof(float)*n);
	cudaMemcpyToSymbol(coord_y,&lig.coord_y[0],sizeof(float)*n);
	cudaMemcpyToSymbol(coord_z,&lig.coord_z[0],sizeof(float)*n);
	cudaMemcpyToSymbol(elec_q, &lig.elecq[0]  ,sizeof(float)*n);
	cudaMemcpyToSymbol(type,   &lig.type[0]   ,sizeof(float)*n);

	t=lig.torCount;
	cudaMemcpyToSymbol(torCount,&t,sizeof(int));

	int axi[t][2];
	unsigned long moving[n];
	bitset<MAX_ATM> bitTor[n];
	bool init_ignore[n];

	float h_axn_x[t];
	float h_axn_y[t];
	float h_axn_z[t];


	for (int i=0;i<n;i++)
	{
		bitTor[i].reset();
		moving[i]=0;
		if (lig.ignore.test(i))
			init_ignore[i]=true;
		else init_ignore[i]=false;
	}

	cudaMemcpyToSymbol(ignore, init_ignore ,sizeof(bool)*n);

	
	map < pair<int, int>, bitset<MAX_ATM> >:: iterator it;


	for (int p=0;p<lig.torCount;p++)
	{
		axi[p][0]=lig.base[p].first;
		axi[p][1]=lig.base[p].second;

		for (int i=0;i<n;i++)
		{
			if (lig.torsion[p].test(i))
				bitTor[i].set(p);
		}

		h_axn_x[p]=(lig.coord_x[axi[p][0]]-lig.coord_x[axi[p][1]]);
		h_axn_y[p]=(lig.coord_y[axi[p][0]]-lig.coord_y[axi[p][1]]);
		h_axn_z[p]=(lig.coord_z[axi[p][0]]-lig.coord_z[axi[p][1]]);
		float dist=sqrt(h_axn_x[p]*h_axn_x[p]+h_axn_y[p]*h_axn_y[p]+h_axn_z[p]*h_axn_z[p]);
		h_axn_x[p]/=dist;
		h_axn_y[p]/=dist;
		h_axn_z[p]/=dist;

	}

	for (int i=0;i<n;i++)
	{
		moving[i]=bitTor[i].to_ulong();
	}

	cudaMemcpyToSymbol(torsion,axi,sizeof(float)*2*t);
	cudaMemcpyToSymbol(move,moving,sizeof(unsigned long)*n);
	cudaMemcpyToSymbol(axn_x,h_axn_x,sizeof(float)*t);
	cudaMemcpyToSymbol(axn_y,h_axn_y,sizeof(float)*t);
	cudaMemcpyToSymbol(axn_z,h_axn_z,sizeof(float)*t);
}

__device__ void quat_rot(float x, float y, float z, float w, float *r)
{
	float tx=x+x;
	float ty=y+y;
	float tz=z+z;
	float twx=w*tx;
	float omtxx=1.0-x*tx;
	float txy=y*tx;
	float txz=z*tx;
	float twy=w*ty;
	float tyy=y*ty;
	float tyz=z*ty;
	float twz=w*tz;
	float tzz=z*tz;

	r[0]=1.0-tyy-tzz;
	r[1]=txy-twz;
	r[2]=txz+twy;
	r[3]=txy+twz;
	r[4]=omtxx-tzz;
	r[5]=tyz-twx;
	r[6]=txz-twy;
	r[7]=tyz+twx;
	r[8]=omtxx-tyy;
}

/*
__global__ void readLig()
{
	int id=threadIdx.x;
	if (id==0) printf("%d atoms in total\n",atomCount);

	printf("thread %d\t[%7.4f %7.4f %7.4f]\t@%7.4f\t%x\n",id,coord_x[id],coord_y[id],coord_z[id],elec_q[id],move[id]);

//
	for (int i=torCount-1;i>=0;i--)
	{
		if (move[id]&(1<<i)!=0) printf("1");
		else printf("0");
	}
//
}
*/

__device__ void transform(float *out, float x, float y, float z, float * fr)
{
	// index is the index of the atom
	// (x,y,z) is the translation vector
	// fr is the rotational factor vector
	// out(x,y,z) is the coordinates of output

	float cdx=out[0],cdy=out[1],cdz=out[2];

	out[0]=cdx*fr[0]+cdy*fr[3]+cdz*fr[6]+x;
	out[1]=cdx*fr[1]+cdy*fr[4]+cdz*fr[7]+y;
	out[2]=cdx*fr[2]+cdy*fr[5]+cdz*fr[8]+z;

	return;
}

__device__ void twist(float *out, int id, float *t)
// make sure twist is done before any transform
// twist is on the original coordinates
{
	for (int i=0;i<torCount;i++)
	{
//		printf("in twist(), move[%d]=%d\n",id, move[id]);
		if (((move[id]&(1<<i))!=0)&&(t[i]!=0))
		{
			float theta_half=t[i];
			float a=sin(theta_half)*axn_x[i];
			float b=sin(theta_half)*axn_y[i];
			float c=sin(theta_half)*axn_z[i];
			float w=cos(theta_half);
			float fr[9];
			quat_rot(a,b,c,w,fr);
/*
			if (id==29)
			{
				printf("in twist(), %f %f %f %f\n",a,b,c,w);
			}
*/
			int ori=torsion[i][0];
			float ori_x=coord_x[ori];
			float ori_y=coord_y[ori];
			float ori_z=coord_z[ori];
			float cdx[3];
			cdx[0]=out[0]-ori_x;
			cdx[1]=out[1]-ori_y;
			cdx[2]=out[2]-ori_z;
			transform(cdx,ori_x,ori_y,ori_z,fr);
			out[0]=cdx[0];
			out[1]=cdx[1];
			out[2]=cdx[2];
		}
	}
}

__device__ float eIntCalc(int index_dist,int atom1, int atom2)
{
	float local=0;
	local+=(elec_q[atom1]*elec_q[atom2]*tex1Dfetch(r_eps_ruler,index_dist+0.5f));
	if (index_dist<=800)
	{
		local+=(tex3D(vdw_tables,index_dist+0.5f,type[atom1]+0.5f,type[atom2]+0.5f));
		float desol=VOL[type[atom1]]*(SOLPAR[type[atom2]]+0.01097f*fabs(elec_q[atom2]))+VOL[type[atom2]]*(SOLPAR[type[atom1]]+0.01097f*fabs(elec_q[atom1]));
		local+=(tex1Dfetch(sol_ruler,index_dist+0.5f)*desol);
	}
	return local;
}

__global__ void trySample(short *nonbonded)
{
	__shared__ float now_x[MAX_ATM];
	__shared__ float now_y[MAX_ATM];
	__shared__ float now_z[MAX_ATM];
	__shared__ float energy[MAX_ATM];

	int id=threadIdx.x;

	float pos[3]={coord_x[id],coord_y[id],coord_z[id]};

	//get torsions by blockIdx

	float tor_angle[]={ 1.460f, -0.385f,  0.895f, -1.422f,  1.258f,  0.519f};

	//apply torsions to both ligand and flex

	twist(pos,id,tor_angle);
	now_x[id]=pos[0];
	now_y[id]=pos[1];
	now_z[id]=pos[2];

	__syncthreads();

	//calculate internal energy
	energy[id]=0;

		for (int i=id+1;i<(id<ligandCount?ligandCount:atomCount);i++)
		{
			if (nonbonded[id*atomCount+i])
			{
				
				int index_dist=(int)(100*sqrt((now_x[id]-now_x[i])*(now_x[id]-now_x[i])+
					(now_y[id]-now_y[i])*(now_y[id]-now_y[i])+
					(now_z[id]-now_z[i])*(now_z[id]-now_z[i])));
				/*
				energy[id]+=(elec_q[id]*elec_q[i]*tex1Dfetch(r_eps_ruler,index_dist+0.5f));
				if (index_dist<=800)
				{	
					energy[id]+=(tex3D(vdw_tables,index_dist+0.5f,type[i]+0.5f,type[id]+0.5f));
					float desol=VOL[type[id]]*(SOLPAR[type[i]]+0.01097f*fabs(elec_q[i]))+VOL[type[i]]*(SOLPAR[type[id]]+0.01097f*fabs(elec_q[id]));
					energy[id]+=(tex1Dfetch(sol_ruler,index_dist+0.5f)*desol);
				}*/
				energy[id]=eIntCalc(index_dist,id,i);
			}
		}

//	printf("thread %3d: %10.4f\n",id,energy[id]);
	float internal=0;
	__syncthreads();
	if (id==0)
	{
		for (int i=0;i<atomCount;i++)
			internal+=energy[i];
		printf("total %10.4f\n",internal);
	//	printf("after torsion\n");
	}
	
	//spawn for translation and oritation // oritation may be pre-generated, http://healpix.jpl.nasa.gov/index.shtml

	//for testing, fix oritation
	float fr[9];
	quat_rot(0.529539f, -0.490877f, -0.339952f, 0.602545f,fr);
	float x=2.543500f, y=6.025500f, z=-7.966300f;
	float bestScore=10000;
	float bestTO[]={0,0,0};
	float hold_x=now_x[id],hold_y=now_y[id],hold_z=now_z[id];
	float intra=energy[id];
	float eq=elec_q[id],abeq=abs(eq);

for (int tx=0;tx<100;tx++)
for (int ty=0;ty<100;ty++)
for (int tz=0;tz<100;tz++)
{
	x=2.0f+((float)tx)*0.01f;
	y=5.5f+((float)ty)*0.01f;
	z=-8.5f+((float)tz)*0.01f;
	//apply translation and oritation to ONLY ligand

	energy[id]=intra;
	pos[0]=hold_x;
	pos[1]=hold_y;
	pos[2]=hold_z;

//	printf("{%.3f %.3f %.3f} thread %2d [%7.3f %7.3f %7.3f]\n",x,y,z,id,pos[0],pos[1],pos[2]);
	__syncthreads();

	if (id<(ligandCount))
		transform(pos,x,y,z,fr);
	now_x[id]=pos[0];
	now_y[id]=pos[1];
	now_z[id]=pos[2];
	
//	printf("thread %2d [ %7.3f %7.3f %7.3f ]\n",id,now_x[id],now_y[id],now_z[id]);

	__syncthreads();

	//calculate the inter-molecule energy
	
	float vdw=0,dsol=0,non=0,ele=0;

	float lkid_x=(pos[0]-min_coord[0])/spacing;
	float lkid_y=(pos[1]-min_coord[1])/spacing;
	float lkid_z=(pos[2]-min_coord[2])/spacing;

	arrayTexFetch(type[id],lkid_x+0.5f,lkid_y+0.5f,lkid_z+0.5f,vdw);
	arrayTexFetch(IDX_D,lkid_x+0.5f,lkid_y+0.5f,lkid_z+0.5f,dsol);
	arrayTexFetch(IDX_E,lkid_x+0.5f,lkid_y+0.5f,lkid_z+0.5f,ele);

	if (ignore[id])
	{
		non=0.0;
		ele=0.0;
	}
	else
	{
//		non=vdw+dsol*abs(elec_q[id]);
		non=vdw+dsol*abeq;
//		ele*=elec_q[id];
		ele*=eq;
	}

	energy[id]+=non+ele;
	
	if (id<ligandCount)
	{
		for (int i=ligandCount;i<atomCount;i++)
		{
			int index_dist=(int)(100*sqrt((now_x[id]-now_x[i])*(now_x[id]-now_x[i])+
				(now_y[id]-now_y[i])*(now_y[id]-now_y[i])+
				(now_z[id]-now_z[i])*(now_z[id]-now_z[i])));
			/*
			energy[id]+=(elec_q[id]*elec_q[i]*tex1Dfetch(r_eps_ruler,index_dist+0.5f));
			if (index_dist<=800)
			{	
				energy[id]+=(tex3D(vdw_tables,index_dist+0.5f,type[i]+0.5f,type[id]+0.5f));
				float desol=VOL[type[id]]*(SOLPAR[type[i]]+0.01097f*fabs(elec_q[i]))+VOL[type[i]]*(SOLPAR[type[id]]+0.01097f*fabs(elec_q[id]));
				energy[id]+=(tex1Dfetch(sol_ruler,index_dist+0.5f)*desol);
			}
			*/
			energy[id]=eIntCalc(index_dist,id,i);
		}
	}
	__syncthreads();

//	printf("thread %2d [ %7.3f %7.3f %7.3f ] %10.3f %10.3f %10.3f %10.3f %10.3f\n",id,now_x[id],now_y[id],now_z[id],vdw,ele,dsol,eq,abeq);
	
	float total=0;
	if (id==0)
	{
		for (int i=0;i<atomCount;i++)
			total+=energy[i];
//		printf("total %10.4f @ [ %7.3f %7.3f %7.3f ]\n",total,x,y,z);

		if (total<bestScore)
		{
		//	printf("better %10.4f @ [ %7.3f %7.3f %7.3f ]\n",total,x,y,z);
			bestScore=total;
			bestTO[0]=x;
			bestTO[1]=y;
			bestTO[2]=z;
		}
	}
	__syncthreads();
}
	if(id==0)
	{
		printf("final best %10.4f @ [ %7.3f %7.3f %7.3f ]\n",bestScore,bestTO[0],bestTO[1],bestTO[2]);
	}

}
/*
__global__ void tryMove(float x, float y, float z, float a, float b, float c, float w, float *t, short *nonbonded)
{
	__shared__ float now_x[MAX_ATM];
	__shared__ float now_y[MAX_ATM];
	__shared__ float now_z[MAX_ATM];

	__shared__ float energy[MAX_ATM];

	int id=threadIdx.x;
	
	float out[3]={coord_x[id],coord_y[id],coord_z[id]};
	float fr[9];
	float theta_half,ta,tb,tc,tw;
	int ori;
	float ori_x,ori_y,ori_z;
	float tcd[3];

//	printf("%f,%f,%f,%f,%f,%f\n",t[0],t[1],t[2],t[3],t[4],t[5]);

//	printf("thread %d\t%lx\n",id,move[id]);

	twist(out,id,t);
//	
	for (int i=0;i<torCount;i++)
	{
		printf("out twist(), move[%d]=%d\n",id,)
		if (((move[id]&(1<<i))!=0)&&(t[i]!=0))
		{
			theta_half=t[i];
			ta=sin(theta_half)*axn_x[i];
			tb=sin(theta_half)*axn_y[i];
			tc=sin(theta_half)*axn_z[i];
			tw=cos(theta_half);
			quat_rot(ta,tb,tc,tw,fr);

			if (id==29)
			{
				printf("in twist(), %f %f %f %f\n",ta,tb,tc,tw);
			}

			ori=torsion[i][0];
			ori_x=coord_x[ori];
			ori_y=coord_y[ori];
			ori_z=coord_z[ori];
			tcd[0]=out[0]-ori_x;
			tcd[1]=out[1]-ori_y;
			tcd[2]=out[2]-ori_z;
			transform(tcd,ori_x,ori_y,ori_z,fr);
			out[0]=tcd[0];
			out[1]=tcd[1];
			out[2]=tcd[2];
		}
	}
//
	quat_rot(a,b,c,w,fr);

	if (id<(ligandCount))
		transform(out,x,y,z,fr);

	now_x[id]=out[0];
	now_y[id]=out[1];
	now_z[id]=out[2];

	__syncthreads();
	
	float vdw=0,dsol=0,non=0,ele=0;

	float lkid_x=(out[0]-min_coord[0])/spacing;
	float lkid_y=(out[1]-min_coord[1])/spacing;
	float lkid_z=(out[2]-min_coord[2])/spacing;

	arrayTexFetch(type[id],lkid_x+0.5f,lkid_y+0.5f,lkid_z+0.5f,vdw);
	arrayTexFetch(IDX_D,lkid_x+0.5f,lkid_y+0.5f,lkid_z+0.5f,dsol);
	arrayTexFetch(IDX_E,lkid_x+0.5f,lkid_y+0.5f,lkid_z+0.5f,ele);

	if (ignore[id])
	{
		non=0.0;
		ele=0.0;
	}
	else
	{
		non=vdw+dsol*abs(elec_q[id]);
		ele*=elec_q[id];
	}

	float e_internal=0,e_lig_int=0,e_flex_int=0,e_lig_flex=0;

	for (int i=id+1;i<atomCount;i++)
	{
//		if ((id<ligandCount&&i>=ligandCount) || (id>=ligandCount&&i<ligandCount))
//		{
//			printf("%d - %d skipped\n",id,i);
//			continue;
//		}
		if (nonbonded[id*atomCount+i])
		{
			float dist=sqrt((now_x[id]-now_x[i])*(now_x[id]-now_x[i])+
				(now_y[id]-now_y[i])*(now_y[id]-now_y[i])+
				(now_z[id]-now_z[i])*(now_z[id]-now_z[i]));

//			if (dist>=8.0f)	continue;

//			float t1,t2,t3;
			int index_dist=(int)100*dist;
			float desol=0;

			//elec
//			float q1q2=elec_q[id]*elec_q[i];
			e_internal+=(elec_q[id]*elec_q[i]*tex1Dfetch(r_eps_ruler,index_dist+0.5f));	//the 48.647321795=332*0.1465;
			//vdw
			if (dist<8.0f)
			{
			e_internal+=(tex3D(vdw_tables,index_dist+0.5f,type[i]+0.5f,type[id]+0.5f));
			//desol:
			desol=VOL[type[id]]*(SOLPAR[type[i]]+0.01097f*fabs(elec_q[i]))+
						VOL[type[i]]*(SOLPAR[type[id]]+0.01097f*fabs(elec_q[id]));
			e_internal+=(tex1Dfetch(sol_ruler,index_dist+0.5f)*desol);
			}

//			printf("atom %2d [%d] with atom %2d [%d] |%d| distance:%11.3f(%4d) elec:%8.3f vdw:%8.3f desol:%8.3f\n",id,type[id],i,type[i],nonbonded[id*atomCount+i],dist,index_dist,t1,t2,t3);

			float temp=elec_q[id]*elec_q[i]*tex1Dfetch(r_eps_ruler,index_dist+0.5f);
			if (dist<8.0f) temp+=tex3D(vdw_tables,index_dist+0.5f,type[i]+0.5f,type[id]+0.5f)+tex1Dfetch(sol_ruler,index_dist+0.5f)*desol;

			if ((id<ligandCount&&i>=ligandCount) || (id>=ligandCount&&i<ligandCount))
			{
				e_lig_flex+=temp;
			}
			else if (id<ligandCount&&i<ligandCount)
			{
				e_lig_int+=temp;
			}
			else if (id>=ligandCount&&i>=ligandCount)
			{
				e_flex_int+=temp;
			}
		}
	}

	energy[id]=non+ele+e_internal;

	__syncthreads();
	
//	printf("thread %d\t%d[%7.3f %7.3f %7.3f]@%d\t%7.2f ,%7.2f\n",id,ignore[id],out[0],out[1],out[2],type[id],non,ele);
	printf("thread %d\t[%7.3f %7.3f %7.3f ]@\t%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f\n",id,out[0],out[1],out[2],non,ele,e_internal,e_lig_int,e_flex_int,e_lig_flex);

	if (id==0)
	{
		float total=0;
		for (int i=0;i<atomCount;i++)
		{
			total+=energy[i];
		}
		printf("total %7.3f\n",total);
	}
}
*/
extern "C" void launch_kernel(int count, float x, float y, float z, float a, float b, float c, float w, float *theta_half,short * nonbonded)
{
//	tryMove <<< 1, count>>>(x,y,z,a,b,c,w,theta_half,nonbonded);
	trySample <<< 1, count >>> (nonbonded);
}
/*
__global__ void test_k(float * out,float x, float y, float z)
{
//	out[0]=tex3D(tex,x+0.5f,y+0.5f,z+0.5f);
	float x0=(float)((int)x),y0=(float)((int)y),z0=(float)((int)z);
	float x1=x0+1, y1=y0+1,z1=z0+1;
	float dx=(x-x0),dy=(y-y0),dz=(z-z0);
	float vdw=0,dsol=0;

	out[0]=tex3D(tex_C,x0+0.5f,y0+0.5f,z0+0.5f);
	out[1]=tex3D(tex_C,x0+0.5f,y0+0.5f,z1+0.5f);
	out[2]=tex3D(tex_C,x0+0.5f,y1+0.5f,z0+0.5f);
	out[3]=tex3D(tex_C,x0+0.5f,y1+0.5f,z1+0.5f);
	out[4]=tex3D(tex_C,x1+0.5f,y0+0.5f,z0+0.5f);
	out[5]=tex3D(tex_C,x1+0.5f,y0+0.5f,z1+0.5f);
	out[6]=tex3D(tex_C,x1+0.5f,y1+0.5f,z0+0.5f);
	out[7]=tex3D(tex_C,x1+0.5f,y1+0.5f,z1+0.5f);

	float i1=out[0]*(1-dz)+out[1]*dz;
	float i2=out[2]*(1-dz)+out[3]*dz;
	float j1=out[4]*(1-dz)+out[5]*dz;
	float j2=out[6]*(1-dz)+out[7]*dz;

	float w1=i1*(1-dy)+i2*dy;
	float w2=j1*(1-dy)+j2*dy;

	out[8]=w1*(1-dx)+w2*dx;

	arrayTexFetch(5  ,x+0.5f,y+0.5f,z+0.5f,vdw);
	arrayTexFetch(300,x+0.5f,y+0.5f,z+0.5f,dsol);

	out[8]=vdw+dsol*0.096;

//	out[9]=tex3D(tex,x+0.5f,y+0.5f,z+0.5f);
	arrayTexFetch(200,x+0.5f,y+0.5f,z+0.5f,out[9])
}
*/
/*
extern "C" void launch_kernel(int count)
{
	readLig <<< 1, count >>> ();
}
*/


/*
__global__ void test_internal(short *nonbonded)
{
	__shared__ float internal_energy[MAX_ATM];
	int id=threadIdx.x;
//	if (id==0) printf("%d %d %d\n",ligandCount,flexCount,atomCount);
//	if (id!=0) return;
	float e_internal=0;

//	for (int i=0;i<999;++i)
//	{
//		printf("vdw_table[%d][%d][%4d]\t%10.3f\n",type[id],5,i,tex3D(vdw_tables,i+0.5f,5+0.5f,type[id]+0.5f));
//		printf("r_eps_ruler[%4d]\t %10.3f\n",i,tex1Dfetch(r_eps_ruler,i+0.5f));
//		printf("sol_ruler[%4d]\t %10.3f\n",i,tex1Dfetch(sol_ruler,i+0.5f));
//	}


	for (int i=id+1;i<atomCount;++i)
	{
		if ((id<ligandCount&&i>=ligandCount) || (id>=ligandCount&&i<ligandCount))
		{
//			printf("%d - %d skipped\n",id,i);
			continue;
		}
		if (nonbonded[id*atomCount+i])
		{
			float dist=sqrt((coord_x[id]-coord_x[i])*(coord_x[id]-coord_x[i])+
				(coord_y[id]-coord_y[i])*(coord_y[id]-coord_y[i])+
				(coord_z[id]-coord_z[i])*(coord_z[id]-coord_z[i]));

//			if (dist>=8.0f)	continue;

			float t1,t2,t3;
			int index_dist=(int)100*dist;

			//elec
			float q1q2=elec_q[id]*elec_q[i];
			e_internal+=(t1=elec_q[id]*elec_q[i]*tex1Dfetch(r_eps_ruler,index_dist+0.5f));	//the 48.647321795=332*0.1465;
			//vdw
			e_internal+=(t2=tex3D(vdw_tables,index_dist+0.5f,type[i]+0.5f,type[id]+0.5f));
			//desol:
			float desol=VOL[type[id]]*(SOLPAR[type[i]]+0.01097f*fabs(elec_q[i]))+
						VOL[type[i]]*(SOLPAR[type[id]]+0.01097f*fabs(elec_q[id]));
			e_internal+=(t3=tex1Dfetch(sol_ruler,index_dist+0.5f)*desol);

			printf("atom %2d [%d] with atom %2d [%d] |%d| distance:%11.3f(%4d) elec:%8.3f vdw:%8.3f desol:%8.3f\n",id,type[id],i,type[i],nonbonded[id*atomCount+i],dist,index_dist,t1,t2,t3);
		}
	}
	internal_energy[id]=e_internal;

//	printf("thread %3d\t%.3f\n",id,internal_energy[id]);

}
*/
/*
extern "C" void launch_kernel(int count, short * nonbonded)
{
	test_internal <<< 1, count >>> (nonbonded);
}
*/
/*
extern "C" void launch_kernel(float * d_out, float x, float y, float z)
{
	test_k <<<1,1>>> (d_out,x,y,z);
}
*/

